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ABSTRACT 

We describe cosmological galaxy formation simulations with the adaptive mesh refinement code 
Enzo that incorporate a star formation prescription regulated by the local abundance of molecular 
hydrogen. We show that this H2 -regulated prescription leads to a suppression of star formation in 
low mass halos (Mh < 10^^ M©) at z > 4, alleviating some of the dwarf galaxy problems faced by 
theoretical galaxy formation models. H2 regulation modifies the efficiency of star formation of cold 
gas directly, rather than indirectly reducing the cold gas content with "supernova feedback". We 
determine the local H2 abundance in our most refined grid cells (76 proper parsec in size at 2: = 4) 
by applying the model of Krumholz, McKee, & Tumlinson, which is based on idealized ID radiative 
transfer calculations of H2 formation-dissociation balance in ~ 100 pc atomic-molecular complexes. 
Our Il2-regulated simulations are able to reproduce the empirical (albeit lower z) Kennicutt-Schmidt 
relation, including the low Egas cutoff due to the transition from atomic to molecular phase and the 
metallicity dependence thereof, without the use of an explicit density threshold in our star formation 
prescription. We compare the evolution of the luminosity function, stellar mass density, and star 
formation rate density from our simulations to recent observational determinations of the same at 
z = 4 — 8 and find reasonable agreement between the two. 

Subject headings: cosmology: theory - galaxies: dwarfs - galaxies: formation - galaxies: halos - 
methods: numerical 



1. INTRODUCTION 

The A CDM paradigm of c o smological structure for- 
mation ( |White fc Rees||1978[ [Blumenthal et al.||1984 ) 
has been tremendously successful at explaining the large 
scale, statistic al features of the dist ribution of matter in 
our universe (Springel et al. 2006). At the same time, 
it is also very clear that the mapping from dark mat- 
ter halos to their baryonic components, to the proper- 
ties of galaxies embedded within the halos, is far from 
straightforward and currently poorly understood. One 
prominent example of this lack of understanding is the 
fact that the cosmic mass-to-light relation is neither con- 
stant nor monotonic, and instead exhibits a minimum 
at a galaxy mass of 10^^ Mq (e.g. Conroy & Wechsler 
|2009). Evidently some unknown processes are inhibiting 
efficient star formation on both larger and smaller mass 
scales. 

The focus of this paper is on the low mass end, the 
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dwarf galaxies. There are (at least) two dwarf galaxy 
problems, which may or may not have the same ex- 
planation. The first of these is the well-known "Miss- 
ing Satellites Problem" 



Kauffmann et al. 1993 Klypin 
et al.||199^ [Moore et^.||19 99), which refers to the dis- 
crepancy between the relatively small number of satel- 



crepancy between tne relatively 
lite galaxies known to be orbiting the Milky Way and 
M31 (~ 20 around each) and the vastly larger number 
of dark matter subhalo satellites predicted from dark- 
matter-only cosmological numerical simulations (> 10^ 
in the latest simulations, [ Diemand et al.| 2008[ jStadel 
|et al. 1 ^2009; Springel et al. 2008j. Photo ' " 



leating from 



the meta-galactic UV background will prevent gas from 



sive subhalos (lEfstathiou 1992 



(lEfsta t: 



all but the most 


mas- 


Kauffmann et al. 


1993 


so the number oJ 


[ dark 



Bullock et al. 

matter subhalos that should be able to host a luminous 
component, either because they collapsed prior to reion- 
ization or because they reached a sufficiently large mass 
thereafter, exceeds the current census of dwarf satellite 



galaxies by at least one order of magnitude ( Madau et al 
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2008). The interpretation of this discrepancy is further 
comphcated by interactio ns with the host galaxy. Ram 
ressure stripping (jMayer et al.|2006), as weh as resonant 

interactions ( [Gnedin 
halo or stellar 



jD'Ongh 
|et al. 1999 

aisk could modify the abundance and properties of galac- 
tic satellite galaxies. 

The second dwarf galaxy problem occurs in the field, 
and is exemplified by the apparent inability of virtually 
all theoretical models of galaxy formation to match the 
abundance of low stellar mass galaxies at z > 0. Semi- 
analytic galaxy formation models (SAMs), for example, 
are able to match the observed stellar mass function in 
the local universe {z = 0) by judiciously tuning their 
AGN, photo- ionization, and supernova feedback param- 
eters, but these same models predict an abundance of 
< 10^^ Mq galaxies at higher redshifts that exceeds 
t he observational constraints by an order of magnitude 



(iFontanot e t al.l 
et al.||2010j ). H; 



2009 



2009 Cirasuolo 



tions race similar 



Ostriker|200"6] 

problem is tin 



Marchesini et al 

ydro dyna mic galaxy tormatTon Simula- 
problems ([Nagamine et al.|2006 
Loi fc Nagamine|201l| ). A closely 



Choi 



Gen fc 

related 

problem is the inability ot numerical simulations and 
SAMs to match the low values of the stellar mass fraction 
and its strongly decreasing trend with halo mass, as in- 
ferred from observations and semi-emp irical approaches 
(|Guo et al ||2010[ |Avila-Reese et al. | [2Ql l ). 

Besides U V photo-heating, stellar feedback, in the form 
of ene rgy or m omentum injection from supernovae explo- 
sions (iDekel fc Silk 1986; Efstathiou„2000), stehar winds 



adia 



Norman fc Silk 1980; McKee 1989) , or radiation pressure 
Xrumholz fc Matzner,2009j ,Murray et al.l2010l [Hopkins 

Dfain till 



et al.||201l[ ), is commonly invoked to explain the reduc 
tion m star fo rmation efficien cy in low ma ss halo s. In 
SAMs (iGole et al. 2000; Benson et al.|[2QQ 3; Some rville 
et al . 2008; Wan g eriTj | 2008| ) this~type of feedbacFTs 
3cr 



typically modeled as a reduction of the cold gas reservoir 
available for star formation, with an efficiency propor- 
tional to some power law of the galactic disk circular 
velocity. 

In direct cosmological numerical simulations, stellar 
feedback remains subgrid physics even for today's state- 
of-the-art computational efforts. Gosmological zoom-in 
simulations of individual galaxies h ave reached tens of 
parsec resolution (Gnedin et al. 200 ^ [Geverino fc Klypin 



2009^ 'Governato et al. 2010; Agertz et al. 2 011[ II^aucKer 



t^iguere fc Keres 2011; Guedes et al. 2011J , but full-box 



behind dOcvirk et al. |2008 Schaye et al.||2010 


lOppen- 


heimer et al. 


mOl IGhoi fc iNagaminel 1:^0111 


I^aucher- 


Giguere et al. 


20 11 p. Neither approach is able to resolve 



erates in nature. Instead simulators have turned to feed- 
back prescriptions that are meant to capture the cumu- 
lative effect of supernovae explosions on scales that are 
computationally accessible. A wide variety of such sub- 
grid physics prescriptions have been implemented, rang- 
ing from a simple injection of thermal energy at t he lo 



cation of newly created star particles (Gen fc Ostriker 



1993), often with radiative cooling artificially turned off 



lor some time to preve nt the newly added energy from 
rapidly radiatin g away ( Thacker fc Gouchman 2000 Stin- 
son et al.|2006 '), to attempts at keeping track ol separate 



|Yepes et ah] 1 1997] |Gnedin| [19981 [Springel fc Hernquist 
!003[ ), to direct kinetic teedback , in which m omentum 



kicks are applied to surrounding gas particles fSpringel 
fc H ernquist 2003; Oppenheimer fc Dave 2006; Schaye^ 
et al.||2010[ [Genel et al.||2010| ), which are subsequently 
te mporarily decoupled from hydrodynamic forces (except 



Schaye et al. 2010) in order to allow them to escape 
the star forming region. 

Although improvements in the stellar feedback treat- 
ment have indeed enabled progress in galaxy formation 
simulations, for example the production of quasi-realistic 
disk galaxies from cosmological initial conditions ([Gov- 



ernato et al.|2010{|Agertz et al.|201l| [Guedes et al |2011[ 
Avila-Reese et al.|2011[|Piontek fc Steinmetz|2011| [Brook 



et al. I2012D, the results depend sens itively on the de - 
tails oithe feedback implementations ([Sales et al.|2010|), 
which themselves are often based on ad-hoc assumptions. 
Furthermore, many problems are not completely solved 
by stellar feedback as it is currently implemented. The 
inability of simulations and SAMs, even those includ- 
ing supernova feedback prescriptions, to match the ob- 
served high redshift stellar mass functio ns and sta r for 



mation rate s is on e example (Girasuolo et al. 2010 Ghoi 
fc Nagami ne|2011 ). Another is the challenge ol suppress- 
ing the stellar mass content of low mass halos while si- 
multaneously matching the observed mass-metall icity re- 
lation in Milky Way dwarf satellites (Font et al. 2011| ). 

Important physical processes associated with star lor- 
mation are not captured by current models, and it is 
time to revisit the dwarf galaxy problems in light of new 
understanding of how SF actually occurs in the ISM of 
local galaxies. One promising direction is an improved 
treatment of the chemistry and thermodynamics of the 
interstellar gas that is actually forming stars. In par- 
ticular, spatially resolved observations of local galaxies 
have revealed that star formation correlates much more 
tightly with th e density of molecular gas than total gas 
density (IWonF fc Blitz 2002 ; jKennicutt et al.|2QQ7|[Leroy 
et al.|2008|[Bigiel et al.||2008| . Even though the primary 



cooling agents are lines ol CO or GII (depending on the 
chemical state of the carbon), molecular hydrogen (H2) 
is expected to be good tracer of star fo rmation, even at 
low metallicities ( Krumholz et al.|[2011 ). 

This motivates a star formation prescription that dif- 
ferentiates between the chemical phases of the gas, in 
which the formation of star particles is tied to the local 
abundance of H2, as opposed to the total gas density, 
as is more commonly done in numerical simulations. In- 
deed some SAMs have begun to explore this direction 



p^ 

del P. Lagos et al.[[2Q10[ Krumholz 



cold and hot phases of the subgrid interstellar medium 



( Fu et al][2010[ 

JJekel 2011), and several numerical simulations inclu ding 
H2 p hysics have bee n published (Robertson fc Kr avtsov[ 
[2Q08[ [Gnedin et al.[|2Q09| [Feldmannet al. 2011).njHtir 
now these simulations focused only on individual galax- 
ies, either in isolated disks or in cosmological zoom-in 
simulations. In this work we investigate for the first time 
the effect of an H2-regulated star formation prescription 
in full-box cosmological simulations, with an eye towards 
the statistical distribution of star formation efficiency in 
dwarf galaxies. 

Following the non-equilibrium H2 chemistry in a re- 
alistic and self-consistent manner, including formation 
on dust grains and the radiative transfer of ionizing 
and dissociating radiation, is complicated and expen- 
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sive to i mplement in nume rical galaxy formation sim- 
ulations ( Gnedin et al.|2QQ9 ). Fortunately, analytical ID 
radiative transfer ca lculations assuniing H2 formation- 



dissociation balanc e ( Krumholz et al.|2QQ8 



Kr umholz 2010 
is determined 



have shown that the 



to a good approximation (Krumholz & 

GnedinpOll ) by the HI column density and metallicity 



2009||McKee 
2_ abundance 



of gas on ~ TOO pc scales. As these scales are directly ac 
cessible to us, we can bypass much of the computational 
difficulty associat ed with proper H2 chemistry by imple- 
menting the Krumholz et al. (2009) results in our cos- 
mological simulations. Even so, our use of cosmological 
adaptive mesh refinement prevents our simulations from 
progressing much past z ~ 4 at an acceptable computa- 
tional expense. Yet it is precisely in the early universe, 
at low but non-zero metallicities, that the metallicity- 
dependence of the HI to H2 transition will be most impor- 
tant, and its effect on star formation greatest. Thanks to 
extensive mult i- wavelength surveys (e.g. the Great Ob- 



servatories Origins Deep Survey, Giavalisco et al. 
and deep follow-up observati ons with the Hu IdK 



2004) 



e and 



Spitzer space telescopes (e.g. Stark et al.""2009; "Bouwens 
|et al. 2011; Gonzalez et al. 2011; Labbe et al. 2010a^7we 
are able to make contact with observational constraints 
on the cosmic stellar mass and star formation density 
even at these high redshifts. 

The main driver of this work, then, is to investigate to 
what degree a proper accounting of the H2 abundance in 
star forming gas can take the role that is traditionally as- 
signed to supernova feedback, namely a reduction of the 
star for mation efficiency in low mass dark matter halos 



(see also IGnedin et al.(|2009[ [Gnedin fc Kravtsov|[2QTQ 



hereafter G09; GKlOp This paper is organized as fol 
lows. In §[2] we describe our numerical approach and the 
details of our H2-regulated star formation prescription. 
We present the results of our work in § [3] - [5j We first 
show that our prescription reproduces many of the obser- 
vational features of the star formation scaling relations 
(§ [3|. We then demonstrate that tying star formation 
to the H2 abundance indeed suppresses star formation 
in low mass halos, thereby alleviating the dwarf galaxy 
problems (§[4|. Finally, we present a direct comparison 
to recent observational determinations of the high red- 
shift evolution of the luminosity function, stellar mass 
density, and star formation rate density (§[5|. We sum- 
marize and conclude in § |6] 

2. SIMULATIONS 

We have conducted cosmological AMR hydrodynamics 
simulations using Enzo v2.(Qto follow galaxy formation 
in the early (2; > 4) universe. The computational do- 
main covers a (12.5 Mpc)^ box with a root grid of 256^ 
grid cells. The dark matter density field is resolved with 
256^ particles of mass 3.1 X lO^M©. The box has mean 
over-density of zero, and no additional density fiuctu- 
ation on the scale o f the simulation box ("DC mode", 
Gnedin et al.||2011|) has been applied. Adaptive mesh 



refinement is allowed to occur throughout the entire do- 
main for a maximum of 7 levels of refinement, resulting in 
a maximum spatial resolution of Axr = 76.3 x 5/(1 + 
proper parsec. Mesh refinement is triggered by a grid 
cell reaching either a dark matter mass equal to 4 times 

^ http://code.google.eom/p/enzo/ 



the mean root grid cell dark matter mass, or a baryonic 
mass equal to 8 x 2~^-^^ times the mean root grid cell 
baryonic mass, where / is the grid level. The negative 
exponent in the baryonic refinement mass threshold im- 
plies a super-Lagrangian refinement criterion that results 
in more aggressive refinement at higher resol ution. The 
simulations are initialized at 2^7 =99 with the |Eisenstein| 



& Hu (1999) transfer function, and cosmologica l parame- 



et al 



ters consist ent with the WMAP 7-year results ( Komatsu 
2011): Q.M = 0.265, = 0.735, Qbh^ = 0.02264, 
Tjul, n = 0.963, and as = 0.801. The parameters of 



our suite of simulations are summarized in Table [H 

The equations of hydrodynamics are solved us- 
ing Enzo's imp lementation of the Piecew ise Parabolic 
Method (PPM, [Colella fc Woodward||1984D , a higher or- 
der accurate Godunov scheme. We utilize the recently 

— r 



added HLLC Riemann solver ( Toro et al. 
a fallback scheme to the more dittusive HL 



1 994| ) with 
solver for 



problematic cells, which greatly aids the simulatio ns' sta- 
bility. En zo employs a dual-energy formalism (Bryan 



et al [1995 ), solving for both the internal gas energy and 



total energy separately, to ensure accurate pressures and 
temperatures in hypersonic fiows. 

Our simulations include radiative cooling from both 
primordial and metal enriched g as, as well as photo- 
heating from an optically thin, uniform meta-galactic UV 
background. The primordial gas cooling rates are calcu- 
lated from the ionization states of hydrogen and helium, 
which are followed with a 6-species (H, H+, He, He+, 
He++, and e~)_ 



et aL| |l 

and pnoto- ionization/ excitation rates. The metal cool 
ing is determined from a 5-dimensional table (indepen 
dent variables: density, temperature, electron fraction, 
metallicity, and redshift) of heating and c ooling rate s 
precomputed with the Cl oudy co de (Ferland et al |1998|), 
as described in detail in Smith et al., ( |2008| ) and ]!Smith 
et al.| (2011). For t he Uv background we us ed tne up- 
dated version of the 'Haardt & Madau' (2001*) UV back- 
ground model that ships with version 07.02.01 of Cloudy. 
This model includes the contributions of both quasars 
and galaxies in a redshift dependent manner, and is in 
reasonable agreement with a more recent calculation o f 
the cosmic UV background fFaucher-Giguere et al. 2009 ). 
As is commonly done in Eulerian hydro dyna mic galaxy 



1997 Anninos 



non-equilibrium ch emical network (l Abel| 
ninos et al. 1997), including collisional 



formation simulations (Machacek et al. 2001; Robert son 
Kravtsov 2008; Agertz et al. 2009^ Ceverino fc Klypin 



2009p , we apply an artificial pressure support to cells 
that have reached the maximum refinement level. This 
is necessary in order to stabilize these cells against artifi- 
cial fragmentation, and is supposed to mimic the pressure 
support from turbulent motions below the simulation's 
resolution. In Enzo this support is implemented by in- 
creasing the internal gas energy up to some multiple of 
the value required to make the cell Jeans stable. We 
have set this factor equal to 10, meaning that the Jeans 
length of the highest resolution cells is artificially set to 
10 ~ 3 times the cell width. 



2.1. Star Formation Prescriptions 

The star formation (hereafter SF) prescriptions we 
have implemented are all variations on the basic Schmidt 
law, whereby the local star formation rate (hereafter 
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Table 1 

Summary of the simulations. 



Name 


^final 


Pgas,SF 


^thresh 




[Zfloor] 


Comment 


KT07 


4.0 


tot 


50cm-3 


— 


— 


Krumholz & Tan| (|2007l SF law 


KT07_low 


6.0 


tot 


5 cm"*^ 


— 


— 


lower SF threshold 


KT07_high 


6.0 


tot 


500cm-3 




— 


higher SF threshold 


KMT09 


4.0 


H2 






-3.0 


[Krumholz et al. (2009): 2-phase equilibrium 


KMT09_L8 


6.0 


H2 






-3.0 


one aaaitional rehnement level (maxlevel=8) 


KMT09_FLW1 


5.0 


H2 




1 


-3.0 


KMT09 with uniform LW 


KMT09_FLW10 


5.0 


H2 




10 


-3.0 


background of 


KMT09_FLW100 


5.0 


H2 




100 


-3.0 


increasing 


KMT09_FLW1000 


5.0 


H2 




1000 


-3.0 


intensity 


KMT09_ZF4.0 


6.0 


H2 






-4.0 


lower Zfloor 


KMT09_ZF2.5 


6.0 


H2 






-2.5 


higher Zfloor 


KMT09_ZF2.0 


6.0 


H2 






-2.0 


even higher Zfloor 


KMT09_ZFzlO 


6.0 


H2 






-3.0 


Zfloor at z = 10 


KMT09_Sob 


5.0 


H2 






-3.0 


Sobolev-like approximation of E 


KMT09_SobL8 


6.0 


H2 






-3.0 


KMT09_Sob with maxlevel=8 



Note. — All simulations have the same box size (12.5 Mpc) and were initialized at 2; = 99 with a WMAP7 cosmology: 
Qm = 0.265, Qa = 0.735, = 0.045, h = 0.71, ag = 0.801, and Us = 0.963. The number of dark matter particles 
is 256^ (mp = 3.64 x 10^ M©), and the root grid dimensions are also 256^. We allow up to 7 levels of adaptive mesh 
refinement, except in KMT09_L8 which has one additional level. Pgas,SF indicates whether the SF prescription is tied 
to the total or H2 gas density, rithresh is the minimum density required for SF to occur (KT07 runs only), Jlw/-^mw 
the intensity of the Lyman- Werner background normalized to the Milky Way's value (KMT09_FLW runs only), and 
[Zfloor] = logio(Zfloor/Z0) is the initial seed metallicity applied at z = 9 in the KMT09 simulations. 



SFR) in a grid cell is proportional to its gas density di- 
vided by a SF time scale, 



PSF 



Pgas 



(1) 



In all of our simulations we set the SF time scale equal 
to the local free-fall time, 



U = tf^ = 



32 G /9gas 



(2) 



Krumholz 
"SF 



fc T an 
etticiency 



an d fix the SF effici e ncy t o e 
cal 



per 



(|2007) 

tree 



0.01, as motivated 
who showed that the lo- 
fall time is low (esF = 



0.01) and approximately constant over 
4 orders of magnitude in density. 

In Enzo's standard routines, SF is allowed to occur at 
every time step and in every grid cell that is not further 
refined. Provided the cell fulfills all conditions for SF, 
a fraction of the cell's gas mass is converted into a star 
particle of mass rrip = e pgas {Ax)^{At/U). For highly re- 
fined cells with small time steps At, this commonly result 
in very large numbers of low mass star particles, which 
can dramatically slow down the simulation's progress. 
Applying a mass threshold, below which a star particle 
is simply not created, is undesirable, since it can lead to a 
significant amount of "unfulfilled" SF, although this can 
be remedi ed with a stochastic SF criterion (e.g. Springel 
& Hernquist 2003, and see below). 

lb overcome these difficulties, we have modified Enzo's 
routines to allow SF to occur only once per root grid time 
step and only in cells at the highest refinement level (here 
/ = 7), but with a star parti cle mass propor tional to the 
root grid time step Ato (see Kravtsov]|2003 ) , i.e. 



Ato 



(3) 



Making rrip proportional to the root grid time step 
{Ato Atr) goes a long way towards overcoming the 



problem of large numbers of low mass star particles. Nev- 
ertheless we also enforce a minimum star particle mass of 
^min = 10^ M0, since even Ato can occasionally become 
very small. Below this mass we implement a stochastic 
SF criterion as follows: if rup < mmin, we form a particle 
of mass equal to mmin if a randomly generated number 
is smaller than {rrip /mmin)- 
We consider two distinct classes of SF prescriptions: 

i) Standard SFR (KT07): The SFR is proportional to 

the total gas density divided by the free-fall time, 

3/2 

resulting in a SFR proportional to pgas • We apply a 
density threshold below which SF is not allowed to 
occur, and vary this threshold between values of 5, 
50, and 500 cm-^. 

ii) H2-regulated SFR (KMT09): The SFR is propor- 
tional to the molecular hydrogen density divided by 
the free-fall time determined from the total gas den- 
sity, resulting in a SFR proportional to Pgas- 
The H2 frac tion , /h^ = PHa /P gas ^ is determined fol- 
lowi ng [Krumholz et al.| ( |2009[ ) (more details in sec- 



tion 12^ )71IH3'~we~consideF^o^ the two-phase equi- 
librium model and a range of models with different 
Lyman- Werner H2-dissociating background intensi- 
ties. No density threshold is applied. 

The resolution of our simulations is not sufficient to re- 
solve the formation sites of the ffist generation of stars, 
the so-called population III. In order to capture the metal 
enrichment resulting from the supernova explosions of 
this primordial stellar population, we instantaneously in- 
troduce a metallicity floor of [Zfloor] = logio (Zfloor /Z©) — 
—3.0 at z = 9, as motivated by recent high resolution 
numerical s imulations of the transition from Pop. Ill to 
Pop. II SF ( [Wise et al.||2012 ). This ensures the presence 
of a minimum amount of metals, which seed subsequent 
star formation and further metal enrichment. We discuss 
the sensitivity of our result s to the time and amplitude 
of this metallicity floor in § 4.3 
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2.2. Molecular Chemistry 
To obtain the molecular hydrogen mass fraction 



PH2/Pgas in a given grid cell, we fohow th e 



model develo ped i n Krumholz et al. 

and IMcKee ^ Krum 



et al. 

mode 



(2008^ 



analytical 



Krumholz 



TTni 



( |2QQ9D , 

is based on a radiative transfer calculation of an 
idealized spherical giant atomic-molecular complex, sub- 
ject to a uniform and isotropic Lyman- Werner (LW) ra- 
diation field. The H2 abundance is calculated assum- 
ing formation-dissociation balance. The solution of this 
problem can conveniently be expressed in three lines: 



/h2 



1 



41 + 0.255' 



ln(l + 0.6x + 0.01x^) 



X = 71 



0.6 



^-16.5/ ^H,o' 



(4) 
(5) 
(6) 



where Tc is the dust optical depth of the cloud, (Jd,-2i is 
the dust cross-section per H nucleus to 1000 A radiation, 
normalized to a value of 10~^^ cm~^, 7^-i6.5 is the rate 
coefficient for H2 formation on dust grains, normalized to 
the M ilky Way value of 10"^^-^ cm^ s"^ (IWolfire et al. 
2008 ) , Gq is ambient UV radiation field intensity. 



normalized to the |Draine| ( 1978 ) value for the Milky Way, 
and nH,o is the volume density of H nuclei in units of 
cm~^. Since both ad and 71 are linearly proportional to 
the dust abundance, their ratio is independent of the gas 
metallicity. 



Krumholz et al. (2009) showed that a further simpli- 
fication to the model can be made if the ISM is as- 
sumed to be in two-phase equilibrium between a cold 
neutral medium (CNM) and and a warm neutral medium 
(WNM) (Wolfire et al. 200 3). The assumption of pres- 
sure balance between these two ISM components forces 
the minimum CNM density to be linearly proportional 
to the intensity of the LW radiation field, with only a 
weak dependence on metallicity: 



31 



l + 3.1(Z/ZsNr 



365 



(7) 



where Zsn is the gas phase metall icity in the solar neigh- 
borhood, and we set Zsn = Z© ( Rodriguez fc Delgado- 
|Inglada'201l]) and Z© = 0.0204. Allowing tor the typical 
lJNM density to be somewhat higher than this minimum 
value, n = ^cnm ^min, we get 



X = 2.3 



0'd,-21 



7^_ 



16.5 



-3.1(Z/Zsn) 



0.365 



Z>CNM 



(8) 



which renders /ho 
intensity. As in 

^CNM = 3 and {ad 



completely independ e nt of the LW 
Krumholz & Gnedin (2011) we set 



2l/'/^-J6.5j 

We have conducted one simulation with this two-phase 
equilibrium assumption (KMT09), and a range of simula- 
tions without it, for which we instead specify a spatially 
uniform LW background intensity equal to 1, 10, 100, and 
1000 times the prese nt day Milky W ay value of 7.5 x 10~^ 
LW photons cm-3 ( |Draine| [T978 | (KMT09_FLW1, 
KMT09_FLW1000). We directly apply these prescrip- 
tions to the highest resolution grid cells in our simula- 



tions, whose size {Ax^ = 54.5 x 7/(1 + z) proper par- 
sec) is comparable to the phy sical dime nsions of giant 
atomic-molecular complexes ( |Blitz||1993l ). The dust op- 
tical depth is given by 



Tc = T^/llH CFd 



^0.067 



Z 

Zsn 



p Axy 

IM0 pC-2 



(9) 



= pAxj is the cell's column density, /in = 
g is the mean mass per H nucleus, and 



where S 
2.3 X 10 

we have set the dust cross section per H nucleus to be 
ad = 10-21 (Z/Zsn) cm2. 

The use of Axj in the calculation of S introduces an 
undesirable explicit resolution dependence in o ur algo- 
rithm, the effect of which we investigate in § |3.4| and 14. 2[ 
A So bolev-like approximation, Ssob = p x (p/^p) (e.g 



G09; Krumholz & Gnedin 2011), would remove the ex- 



plicit aipenHence~orTEe~a1'go^^ on the width of the 
finest grid cells. As discussed in detail in the Appendix, 
we have conducted additional simulations using Ssob in 
the KMT09 prescription, and these resulted in slightly 
lower column densities and reduced SF rates. At densi- 
ties relevant for star formation (n > 5 cm~^), the dif- 
ferences between the two approximations of S are small, 
less than 20% in the median with a scatter of 0.26 dex, 
and well within the range of uncertainty of the parame- 
ters of the KMT09 model. The results presented in this 
paper were obtained with the simpler and computation- 
ally less expensive direct cell-based approximation using 
Sceii pAxr. 

Regardless of how the surface density of atomic- 
molecular complexes is calculated, we should expect 
some degree of residual resolution dependence in our sim- 
ulations, since we don't include the necessary stellar feed- 
back physics (see below) that regulates the structure of 
molecular clouds and provides pressure support against 
further collapse. The artificial pressure support men- 
tioned above, which we apply in order to avoid spurious 
fragmentation, is resolution dependent, and hence higher 
numerical resolution (through additional refinement lev- 
els) will always result in further collapse and higher den- 
sities in our simulations. 

The KMT09 model has recently been tested against 
numerical simulations that self-consistently follow 
the formation and destruction of H2 with a non- 
equilibrium chemical network including time-dependent 
and spatially-inhomoge neous 3D radiative transfer of UV 
and ionizing radiation ( [Krumholz fc Gnedin||20lT ). The 
analytical model agrees extremely well with the numeri- 
cal results whenever the metallicity is around 1 per cent 
solar or greater. This agreement holds for both "fixed 
ISM" simulations in which the metallicity and radiation 
field are kept constant and in cosmological simulations in 
which the metallicity and radiation field are computed 
self-consistently. At metallicities below 1 per cent so- 
lar, the analytical model overestimates fn^. Note that 
for the runs without the two-phase equilibrium assump- 
tion, we follow Krumholz & Gnedin (2011) and apply a 
"clumping factor" of 30 to the H2 formation rate (i.e. 
{o'd,-2i/7l-i6.5) = 1/30) to account for unresolved den- 
sity inhomogeneities below our simulations' resolution 
limit. We do not use a clumping factor for the two-phase 
equilibrium model, because the two-phase model is in ef- 



6 



feet a direet estimate of the proper clumping factor. This 
physicahy-motivated clumping model obviates the need 
for an ad-hoc correction. 

2.3. Feedback 

Although we explore in this work to what degree an 
improved treatment of the ISM's chemical state can re- 
place the need for stellar feedback in regulating star 
formation, some form of feedback is necessary even in 
our simulations, if only to enrich the gas with metals 
that promote the formation of molecular hydrogen. For 
this purpose we employ a very simple feedback mech- 
anism that is meant to simultaneously account for the 
mass, metals, and thermal energy deposited by winds 
from massive stars and core-collapse supernovae. The 
feedback is applied instantaneously at the time of for- 
mation of the star particle, is deposited into the / = 7 
grid cell containing the particle, and consists of the fol- 
lowing three components: (i) a fraction Cm — 0.25 of the 
star particle's mass is returned as gas, Amtot = ^m'^p\ 
(ii) the gas is enriched with a metal yield of F = 0.02, 
Am^ = Amtot (^(1 — ^p) + ^p), where is the metal- 
licity of the star particle and the two terms correspond to 
newly and previously enriched ejected material, respec- 
tively; (iii) the thermal energy of the gas cell is increased 
by a fraction esN = 10~^ of the rest-mass energy of the 
newly formed star particle, A£^ = esN vn^c? . 

This feedback implementation is commonly applied in 
cosmological hydrodynamics simulations, but in fact it is 
known to be insufficiently strong. The problem with this 
instantaneous and localized feedback is that the thermal 
energy is applied to grid cells with very high gas density, 
in which the cooling time is very short. The injected en- 
ergy is thus almost immediately lost to radiative cooling, 
and the feedback "fizzles out" ( |Katz|19 92 ) . We acknowl- 
edge the shortcoming of our current feedback implemen- 
tation, but since we are focusing here on the global ef- 
fects of different star formation prescriptions, we defer 
attempts at improving the treatment of feedback to fu- 
ture work. 

Finally, we caution that nothing in our implementa- 
tion of molecular chemistry obviates the need to rely on 
a subgrid model of star formation. Our model should 
be viewed in the same light as other subgrid models in 
wide use, e.^ . the two-phase model of |Springel fc Hern- 



uist 



(2003) or the blast wave model of Stinson et al. 
2UU6p. The main advantage of our new model is that 



it mcorporates an explicit metallicity-dependence, which 
both observations and theory appear to demand, as we 
discuss in more detail below. 

2.4. Ealo population 

We used the HOP halo finder ( [Eisenstein fc Hut|1998] ) 
to identify gravitationally bound dark matter halos in 
our simulations. For every halo we determined the halo 
center (defined as the location of the highest dark matter 
density), the virial radius and corresponding mass (de- 
fined as t he radius enclosing Ayjr ^ times the background 
density ( [Bryan fc Norman 1998| )), and the amount of 
gas (total. Hi, Hii, H2, Hei, Heii, and metals) contained 
within the halo. 

Fig. [1] shows the cumulative mass function o f the sim- 
ulations at z = 4, compared to the .Sheth fc Tormen 
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Figure 1. The cumulative mass function of halos in our simula- 
tions at z = 4. The thick line is for total mass (KT07 and KMT09 
curves are almost indistinguishable, so we only show KT07), and 
the thin lines are for the stellar mass. The dashed line shows the 
Sheth Tormen| ( ||1999, ) mass function fit. The right ordinate gives 
the total number of objects found in our (12.5 Mpc)"^ simulation 
volume. 



(1999) fit. The agreement is remarkably good down to 
M ^ 2x 10^ M0, corresponding to halos with ~ 500 dark 
matter particles. At even lower mass numerical resolu- 
tion effects lead to an artifical suppression in the number 
of halos. A finer root grid resolution and a more aggres- 
sive dark matter based refinement criterion would extend 
the mass function by another order of magnitude in mass 
dO^Shea et aL||2QQ5| . For our purposes, however, a reso- 
lution limit of ^ 2 X 10^ Mq is sufficient, since halos with 
M < 10^ M0 lie below the UV suppression scale . Fur- 
thermore, as we show in Sec. [4j H2 regulated SF leads to 
a suppression of stellar mass in halos with M < 10^^ M©. 
As a brief preview of this effect, we overplot in Fig. [l] the 
cumulative stellar mass functions: the KMT09 simula- 
tion has far fewer low stellar mass (M^ < 10^ M©) galax- 
ies, which is a result of the suppression of SF in halos 
with M < 10^0 Mo . 

3. THE KENNICUTT-SCHMIDT LAW 

A common test for a new numerical SF implementa- 
tion is to compare the simulations to the observational 
Kennicutt-Schmidt (hereafter KS) relation, the empirical 
power law between the SFR surface density an d total gas 
surface density: Hsfr oc Hg^^g, with n ^1.4 (Kennicutt 
19981 [Daddi et al. "2010! [Genzel et al.][2QTQl ). 

In order to facilitate comparisons to observational 
data, we determine column densities by integrating 
through the density fields along the direction perpen- 
dicular to the stellar disk of the simulated galaxies. For 
^SFR we integrate a SFR density defined as 



PSFR 



E 



T* (Ax) 



(10) 



where Ax is the cell width, is the star particle 
mass, and = 10 Myr is the SF averaging time scale, 
which roughly corresponds to observational SFR esti- 
mated from nebular emission lines (Ha, O III) or FIR 
continuum, but is a factor of 5 -10 shorter than estimates 
based on FUV measurements (Kennicutt^ 1 998^ Feldmann| 
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Figure 2. Surface density of SFR, total gas, H2, and metallicity for four representative massive halos (M^ = 2.5 x lO-*^-*^, 2.0 x lO-*^-*^, 1.6 x 
lO-*^-*^, 7.4 X 10-*^° M0) in the KMT09 simulation at z = 4. The projections are calculated as line integrations perpendicular to the disk plane, 
for a 1000 x 1000 pixel grid covering a 10 x 10 kpc region centered on each galaxy. 



|et al.|2QTT|. The sum in Equation 10 is over all star par- 
ticles with age less than r^. 

We calculate Sgas, ^Hs, ^hi, and Ssfr on a 1000 x 
1000 pixel grid covering alO x 10 kpc region centered on 
each galaxy in a sample of 35 of the most massive galax- 
ies in our simulations, chosen to cover a wide range in 
mean metallicity (see § |3.2[ ). We only considered galax- 
ies in which the adaptive mesh refinement reached the 
maximum level (/ = 7). Visualizations of Ssfr, ^gas^ 
, and a density- weighted projection of metallicity are 



in the KMT09 simulation at z = 4. 
The observational KS rela tion, as reported fo r exam - 



pie by 



Daddi et al. 



( |2QTQ ) 



and 



Genzel et al. 



( [20101 ) 



sity . ^ 

shown in Figure [2] for four representative massive galaxies lor references). A reproduction of t 



has been established from a wide variety of data out 
to z > 3, including spatially resolved local z = spi- 
ral galaxies, infrared- selected starbursting galaxies and 
(U)LIRGs, BzK-selected galaxies at z ~ 1.5, "normal" 
star- forming galaxies at z = 1 — 2.3, and starburst ing 
sub-millim eter galaxies ( SMGs ) at z 1 — 3 (see Ke nni-| 



cutt 



( |1998D , [Daddi et al.] ( [2QTQl ), and lGenzel et al., (,2a 



le observational re- 




S_[Mgpc-2] S„3,[Mgpc-2] S„3,[Mgpc-2] 



Figure 3. Comparison of the observational Kennicutt-Schmidt relation from Daddi et aiy] ( |201 0') (left panel) to a standard SF simulation 
with SF density threshold of rithresh = 50 cm~^ (middle) and an H2-regulated SF simulation without any density threshold (right) at 
z = 4. In the simulations the surface densities have been determined from line integrations perpendicular to the galaxies' disk plane, and 
have been smoothed to a resolution of 200 pc (see text for more detail). 
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Figure 4. Comparison of the relation between total SFR and 
total H2 mass for galaxies in the KMT09 simulation at 2: = 4 

iblue triangles) with the observational results from Daddi et al.| 
2010) (their Fig.l) for lower redshift galaxies. The simulation and 
observations have relations with similar slopes, and the simulated 
galaxies lies in between the "sequence of starbursts" and "sequence 
of disks". Higher SFR for a given H2 mass may be expected for 
lower metallicity systems, which require higher densities to allow 
the transition to H2 to occur. 



|2Q1Q[ ) (their Fig.2) is shown in 



lation from Daddi et al.^ 
the left panel of l^lgure [3 

Our simulated galaxies are dwarf galaxies with stellar 
masses less than 10^^ Mq at high redshift ^ > 4, and 
hence not directly comparable to any of these observa- 
tional galaxy samples. Nevertheless, since SF is a local 
process, it makes sense to directly compare our simulated 
KS relation to the observations, as a test of how well our 



SF prescription is performing. The center and right pan- 
els of Figure |3] show this comparison for the KT07 and 
KMT09 simulations at z = 4. Since the 10 pc intrinsic 
resolution of our surface density maps is f iner than in 
most spatially resolved studies to date (e.g. Bigiel et al. 



2QQ8[ |Bolatto et al.|2Ql l), we downgrade the spatial res- 
olution to 2UU pc by spatially averaging with a 20-cell 
boxcar average. We explore the resolution dependence 
of our KS relations in § |3.3| Only pixels with non-zero 
SsFR are plotted. The solid line represents a sliding av- 
erage of log;LO ^gas in logio ^SFR bius of width 0.5. 

Both the amplitude and the slope of the relation are 
qualitatively in very good agreement with the obser- 
vations. S uch agreement has previously been reported 



(|Kravtsov| [20031 J Gnedin fc Kravtsov| |2QlQ[ [2011 
mann et al.(20llp , and we show it here merely to c 



Feld- 



demon- 

strate that our new SF algorithm is valid and that the use 
of our H2-regulated SF pre scription does no t destroy this 



agreement. We note that Kravtsov (2003) showed that 



even a linear local SF law, with a constant SF time scale 
and a density threshold, can, when spatially averaged 
on < kpc scales, result in a super-linear surface density 
relation in agreement with the empirical KS law. Match- 
ing the observed KS law in cosmological simulations is 
thus not a good way to distinguish between different SF 
implementations. 

Another cross check with observational data is pro- 
vided by the demonstration in Fig. |4] that our simulated 
galaxies exhibit a similar relation between their S FR and 
total H2 co ntent as the galaxies in the study by |Daddi 



et al.| ( |2010[ ). Our simulated relation has the same slope, 
and lies m amplitude between what they refer to as the 
"sequence of starbursts" and the "sequence of disks". 
Note that a somewhat higher SFR at a given H2 mass 
may be expected for lower metallicity systems, for which 
the transition to H2 occurs at higher column densities. 

3.1. Low Sgas SF Threshold 

Both KT07 and KMT09 simulations exhibit a drop off 
in SsFR at Sgas ^ 50 pc~^. This feature is believed 
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Figure 5. Simulated mean Kennicutt-Schmidt relations for 
standard SF with three different density thresholds, nthresh — 
5, 50, 500 cm-3, and the H2-regulated SF. The KT07 and KMT09 
lines are from the 2; = 4 output, the KT07_low and KT07_high ones 
from 2; = 6. 



to correspond to saturation in the atomic hydrogen frac- 
tion, with gas becoming primarily molecular at higher 
surface densities ( Bigiel et al.,2008) . Since Ssfr appears 
to be almost independent ot Uhi and instead correlates 
primarily with , this HI saturation point is reflected 
in a kink towards lower SFR in the KS relation. 

In the KT07 simulation, which doesn't account for the 
atomic to molecular hydrogen transition, this cutoff is 
reproduced by the density threshold imposed on the SF. 
In Figure [5] we show how the simulated KS relation de- 
pends on the value of this threshold. As expected the 
relation extends to lower Sgas for the lower threshold 

^thresh 

of S 



= 5 cm case and steepens at a larger value 
^ 100 M0 pc~^ w hen the threshold is higher, 
500 cm"^ (see also Colm et al.|2010 ). Once the 



gas 

'^thresh ~ 

SF is regulated by H2, however, this turnover arises nat- 
urally: although the curve is noisier due to the smaller 
statistics, the relation appears to naturally steepen at 
Sgas ^ 50 Mq pc~^. This confirms the earlier results 
found by GKIO for a set of cosmological zoom-in galaxy 
formation simulations including non-equilibrium H2 for- 
mation and radiative transfer. 

Note that in the KMT09 KS relation (right panel of 
Fig. [3| there are a few points with very low Sgas yet 
non-zero Ssfr^ and a roughly corresponding number of 
points with high Sgas yet values of Ssfr significantly 
below the KS relation. These points arise from star par- 
ticles having wandered out of the high density cells in 
which they were born into a neighboring pixel with much 
lower Ugas- The degree to which this wandering causes a 
smearing in the KS relation depends on the spatial and 
temporal averaging scales employed in calculating Sgas 
and SsFR, and on the stars' velocity dispersion. 

In local normal disk galaxies stellar velocity dispersions 
are typically only 5-10 km s~^, preventing stars from 



Figure 6. Kennicutt-Schmidt relation for total (red), H2 (blue), 
and HI (magenta) surface density from the KMT09 simulation at 
z = 4. The solid lines indicate the mean relations from the simula- 
tions, the shaded regions cover the central 68% scatter (16*^-84*^ 
percentile), and the two da shed lin es are the best-fit ~ ^SFR 
relations reported by Genzel et al.'('2010) for their z = — 3.5 sam- 
ples of "normal" SF galaxies (lower line) and luminous mergers 
(LIRGs/ULIRGs and SMGs) (upper line). 
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Figure 7. The total gas KS relation for simulated galaxy subsam- 
ples split by their mean gas-phase metallicity, from the KMT09 
simulation at 2; = 4. Solid lines indicate the mean relations from 
the simulations, and the shaded regions cover the central 68% scat- 
ter (16^^ — 84*^ percentile). The lower the metallicity, the higher 
the Sgas that is needed for the transition to fully molecular gas, 
which corresponds to the turn over in the KS relation. 
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Figure 8. Comparison of the H2-to-HI ratio versus total neu- 
tral hydrogen column density (Ehi + ^ h?) between SMC data 
smoothed at 200 pc (from Bolatto et al.|2011) and low metallicity 
gas in the KMT09 simulation at 2; = 4 also smoothed at 200 pc. 
The overlap between simulated and observational data indicates 
that our model does an adequate job of capturing the atomic-to- 
molecular transition in real dwarf galaxies. 



wandering over more than a small fraction of the ~ kpc 
smoothing scale in the ~ 10 Myr over which they produce 
significant ionizing luminosities that may be detected in 
Ha or other nebular emission lines. Typical velocity dis- 
persi ons, however , ma y well be larger at higher redshift 
(see jCresci et al. ( 2009, ) for some empirical evidence for 
this), as expected it they are se t by cosmological accre- 



tion ( [Krumholz &: Burkert|20l"Q ) instead of cold disk dy- 
namics. In fact, we tind that averaged on 1 kpc scales 
our galaxies have a mean ID stellar velocity dispersiorP] 
of 44 km s-^ at z=4. In 5% of ah cells (19 in total) 
the ID stellar velocity dispersion exceeds 100 km s~^. 
At these velocities wandering becomes more important, 
and as a result an increased scatter in the KS relation 
should be expected at high redshift. We have verified 
that the amount of scatter at low Sgas increases when we 
reduce the spatial averaging scale of the surface density 
ma ps (see Fig. ^ or increase the temporal SF timescale 
( Feldmann_fc Gnedin||201l{ see also). 
In l^lgure |5] we compare the KS relations for the to- 
tal, atomic, and molecular gas in the KMT09 simulation. 
The Il2 relation does not exhibit a cutoff at low column 
densities, since we have not imposed any explicit density 
threshold in the H2-regulated SF prescription. The H2 
relation is shallower than the total gas one, and its slope 
is in excellent agreement with the observational determi- 
nation of the slope of the molecular gas KS relation by 
[Genzel et al. (2010). These authors studied a population 
of ''normal" star forming galaxies at z — — 3.5 and a 
population of luminous z ~ and 2; ~ 1 — 3.5 mergers 

^ We only include cells containing more than 10 star particles in 
the average. 



(LIRGs/ULIRGs and SMGs) and found that both had 
equal H2 KS slopes of 1.17 with about a 1 dex higher 
normalization for the luminous merger sample. The two 
best-fit relations from their work are shown as dashed 
lines in Figure [6j and our H2 KS relation lies right in be- 
tween the two. The HI relation instead is much steeper 
and doesn't extend much beyond 70 pc~^. The near 
constant offset at low Ugas between the total gas and 
HI KS relations is due to the presence of a signification 
amount of ionized gas on ~ kpc scales. The overall pic- 
ture matches the results reported by GKIO and is in qual- 
itative agree ment with the empirical findings reported by 
Bigiel et al.| ([2Q08), confirming that the origin of the turn 
over in the total gas KS relation is indeed the transition 
from predominantly atomic to fully molecular gas. 

3.2. Metallicity Dependence 

The KS cutoff occurs at somewhat higher Sgas in both 
KT07 and KMT than in the observational data from 



Bigiel et al. ( 2008^ . This can be understood as a result of 
the lower metallicities in our simulated galaxies. We split 
our 35 galaxy sample into three subsamples of different 
mean gas phase metallicities, (Z/Z©) = 0.055, 0.17, and 
0.32. Each subsample contains at least 10 galaxies. The 
highest metallicity galaxy has Z = 0.56 Z©, considerably 
below the closer to sola r metallicity sani ple of 2: = field 
galaxies analyzed by^igiel et al. ( 2008[ ). The lower the 
metallicity of the gas, the higher the total gas surface 
density that is required in order to provide enough LW 
shielding to allow the transition to fully molecular gas. 
The Sgas scale at which the KS cutoff occurs is thus 
expected to scale inversely with metallicity, and this is 
exactly what Figure [7| shows. The cutoff in the KS re- 
lation shifts to progressively higher Ugas for decreasing 
mean metallicity, occuring at > 100 M© pc~^ for the 
lowest metallicity case. A similar trend was previously 
reported by GKIO. Recent observations of star forma- 
tion in the Small Magellanic Cloud (SMC) show that the 
break in the total gas star formation law is indeed shifted 
to higher surface density by a factor of Zct^/Zsm c ^ 5, 
in precisely the manner that our models predict ( Bolatto] 



et al. 20111). The higher Sgas cutoff in the KS relation 
tor low metallicity systems may also be responsible for 
the observational results that Esfr in damped Lyman-a 
(DLA) sys tems ([^^ Ife & Chen 2006) and Lyman-break 
galaxies (Rafelski et al. 2009, 2011) at z 3 appears to 
lie well below the z = KS relation (see GKIO for a more 
in depth exploration of this possibility). 

The success of our model at capturing the physics gov- 
erning the gas phase structure in dwarf galaxies is fur- 
ther demonstrated by Fig. ^ There we compare the 
H2-fraction (Shs/^hi) as a function of total neutral hy- 
drogen gas column (Shi + ^Hs) foi" the low metallicity 
galaxy sample ((Z) = 0.055 Z©) with observational data 
for the SMC that has recently become available. |Bo-| 
(latto et al.| (1201 1 ) have determined a H2 column density 
map of the SMC at ~ 12 pc resolution by combining 
Spitzer IR measurements and radio (ATCA and Parkers) 
HI data. In Fig. [S] we compared their data smoothed on 
~ 200 pc scale (kindly provided by A. Bolatto) to our 
z = 4 KMT09 simulation data smoothed at the same 
scale. There is good overlap between simulation and ob- 
servational data, but in our simulations the distribution 
of points extends both to higher total gas columns and to 
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Figure 9. Dependence of the total gas (left) and H2 (right) KS relation on the smoothing scale A in the KMT09 simulation at z = 4. 
The solid lines show the mean relations and the shaded regions cover the central 68% scatter (for clarity we only show the scatter for A = 
10 and 1000 pc). The mean relation does not show muc h dependence on averaging scale, but the scatter increases towards smaller A, in 
agreement with observations ( [Schruba et al.|2010[ |Liu et al. 2011) . 



lower H2-fractions than probed by Bolatto et al. (2011) 



3.3. Smoothing Scale Dependence 

In Figure [9] we show the dependence of the total gas 
and H2 KS relations on the spatial scale over which the 
data is smoothed. In addition to our fiducial scale of 
200 pc, we present the relations for up to twenty times 
finer and five times coarser smoothing. This range of 
smoothing scales roughly mimics the variations in the 
angular resolution (i.e. beam size) of the radio obser- 
vations used to establish the observational KS relation, 
ranging from spatially resolved measur ements of nearby 
spira l galaxies with sub -kpc resolution (Kennicutt et al. 
[20071 |Bigiel et al.||2008|) to high redshift observations in 
which a large traction of the galaxy is contained in a sin- 
gle beam (Kennicirttl998). We do not see much evidence 
for a smooth scale dependence, in either total gas or H2 
KS relations. The scatter in the relations, however, in- 
creases for s maller smoothing sca l es, in agreemen t with 
observations ( Schruba et al.|[2'010 Liu et aL]|2011 ). 



3.4. Resolution Dependence 

With only seven levels of adaptive mesh refinement, 
our simulations are unable to resolve the true Jeans 
length of the cold, molecular gas in star forming galax- 
ies. As discussed in § [2j we resort to an artificial min- 
imum pressure support in order to stabilize gas cells at 
the highest refinement level against artificial fragmenta- 
tion. This has the undesirable consequence of making 
the results of our simulations somewhat dependent on 
resolution, since additional levels of refinement will al- 
low gas to collapse further and reach higher densities, 
until the resolution becomes adequate to resolve the true 
Jeans length. Unfortunately, additional refinement lev- 
els come at a computational cost. To run the KMT09_L8 
simulation, a clone of KMT09 with one additional refine- 
ment level (/max = 8), down to z = 6 took about 3 times 
as long as the KMT09 run. Note that virtually every 



/ = 7 KMT09 grid was at least partially further refined 
in KMT09_L8. 

In the left panel of Figure [lO] we show a comparison of 
the distribution functions of number density in the max- 
imally refined grid cells in the KMT09 and KMT09_L8 
simulations at 2: = 6. The additional refinement level has 
allowed gas to collapse to higher densities. The mean 
density at / = 7 is 427 cm"^ in the KMT09 simulation, 
but 1420 cm-^ at / = 8 in KMT09_L8. Higher densities 
will lead to larger a SFR and an increase in Esfr- ^gas 
smoothed on ~ kpc scales, however, will be unaffected, 
since it depends only on the total mass enclosed in a 
given kpc scale column, not the local density. It is not 
surprising, then, that the simulated KS relation (right 
panel of Fig. [lO| shows that the KMT09_L8 KS relation 
has somewhatnigher amplitude than in KMT09. 

It is important to note that, despite its higher reso- 
lution, model KMT09_L8 is not necessarily more real- 
istic than KMT09, because it lacks the physics needed 
to properly model molecular clouds at the increased res- 
olution. In the absence of feedback mechanisms other 
than supernovae, increasing resolution allows the gas 
to collapse to ever-higher density, so that the bulk of 
the molecular gas will always reside near the resolution 
limit. However, this behavior is not realistic. In ob- 
served nearby galaxies, the bulk of the molecular clouds 
exist at densities of a few hundred cm~^ rather than a 
few thousand cm ~^ (i.e. clo ser to the mean in KMT09 
than KMT09_L8; |Blitz|1993[) , and molecular cloud prop- 
erties are independent ot gala ctic environment, st rongly 
suggesting internal regulation ( Bolatto et aH20Q8 ). Pos- 
sible mechanisms to provide this regulation m clude ion- 
ized gas pr essure (Matzner 2002 [KrumholFet aL 20061 
Goldbaum et al. 2011|), proto stellar winds ( [INakamura 
& Li 2007; "Wang et al.''2010l) , and radiatio n pre s sure 



2010 



Fall 



(Krumholz &; Matzner 2009; Murray et al^ ^ 

|et al.n2010| (Hopkins et al.||20lf^ none of which are m- 
eluded in our simulations. Lhus we regard the KS law 
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Figure 10. Left: Probability distribution functions of the proper gas number density in the KMT09 and KMT09_L8 simulations at 2; = 6 
for cells at the maximum refinement level (/max = 7 and 8, respectively). Right: The corresponding total gas KS relation. The additional 
refinement level allows gas in KMT09_L8 to reach higher densities, which in turn results in greater SFR densities, and an increase in the 
normalization of the KS relation. 



we obtain from KMT09 as at least as reliable as the one 
we obtain from KMT09_L8. Moreover, this comparison 
yields an important caution: increasing resolution with- 
out a corresponding increase in physics does not neces- 
sarily yield a better result. 

4. STELLAR MASS FRACTION 

We now turn to the effects of an H2-regulated SF pre- 
scription on the stellar content of the dark matter halos 
in our simulations. In Figure [TT] we show plots of the 
stellar mass fraction f^ = M^/Mh against the total halo 
mass Mh = MoM+Mgas+M^ for simulations with a stan- 
dard (KT07) and Hs-regulated (KMT09_ZFzl(f]) SF pre- 
scription. In the normal SF case (left panel), nalos with 
total masses as low as 10^ Mq have been able to form 
a substantial stellar content, with values of f^ not drop- 
ping much below 1%. This is problematic in view of the 
observational dearth of such dwarf galaxies in the local 
universe. Given the high volume density of M > 10^ Mq 
dark matter halos predicted by ACDM structure forma- 
tion, such a high star formation efficiency would vastly 
overproduce the faint end of the field dwarf galaxy lumi- 
nosity function and the abundance of faint Local Group 
dwarf galaxies. 

In Via Lactea II, a collisionless simulati on of the forma- 
tion of a Milky Way's dark matter halo (|Diemand et al. 
[2008 ) , for example, there are ~ 100 halos with an identi- 
nable z = remnant within the host halo's virial volume 
that had a mass exceeding 10^ Mq at some point in their 
evolution. Yet of the Milky Way dwarf satellite galax- 
ies, probably only the two Magellanic Clouds, Fornax, 
and Leo I (and possibly the progenitor of the Sagittarius 
dwarf) have a stellar mass greater than 10^ M©. The vast 
majority of such dark mat ter halos must thus have SF ef- 
ficiencies well below 1% ( Rashkov et al.|2012|) . The same 
is true of faint field galaxies, not just satellites. Com- 

^ We now focus on KMT09_ZFzlO, because the slightly earlier 
metallicity floor (z = 10, instead of z = 9 in KMT09) results in a 
more gradual suppression of (f^), see § |4.3| 



parisons of the local universe galaxy luminosity function 
from SDSS with the predicted DM field halo mass func- 
tion show that at z < 1 the stellar mass fraction must 
have a maximum of ~ 0.01 — 0.03 aroun d 10^^ M^, de- 
creasing sharply towards lower masses (Zheng et al. 2007| 



[Conroy & Wechsler 20091 [G uoet al. 2010; Moster et 
[2010; Behr oozi et al.| 201Q|. 

T^u Jy 



"I'he most commonly invoked mechanisms to lower the 
stellar mass fraction in these low mass galaxies are (i) 
the prevention of gas cooling due to the meta-galactic 
UV background (Efstathiou 1992; Kauffmann et al. 1993, 
[Bullock et al.|20Q0| ) and (ii) stellar feedback f rom massive 



winds and supernoyae ([ L arson ^ ^1974^ Kauffmann et al 
'2002|. Cnir simulations include a 



1993 



Benson et al. 



model lor the U V background, via a spatially uniform, 
optically thin radiation field that accounts f or the UV 
emissivity from both quasars and galaxie^ ("Haard t fc] 
Madau 2001) . By itself this feedback does not appear to 
be able to suppress star formation in 10^ — 10^ ^ halos 
at 2; > 4 in our simulations. As discussed in § [2.3[ we do 
include a supernova feedback prescription in our simula- 
tion, but in a form known to be insufficiently strong to 
effect the necessary reduction in SF efficiency or to repro- 
duce observed galactic outflows. This is exemplified by 
the high stellar mass fractions (f^ ^ 0.06, a factor 5 — 10 
higher than what is allowed by observational constraints 
at much lower redshift, ^ ^ 1) e ven in our m ost massive 
halos (10^^ - 10^^ M0) (see also Avila-Reese et al. 20lT| . 

It's reasonable to expect a similar, or even greater, 
suppression of f^ for lower mass halos, whose shallower 
potential wells should make it easier for winds and su- 
pernovae to expel their gas 



Font et al. 



(2011) 



However, 

find that in order to simultaneously niatch the lummos 
ity function of Local Group dwarf galaxies and their ob- 
served mass-metallicity relation, the strength of stellar 
feedback must saturate in halos with Fmax ^ 65 km s~^ 
5 X 10^^ M0). Furthermore, a physically realistic 

^ The galactic contribution dominates the UV background at 
this early epoch. 
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Figure 11. Stellar mass fraction vs. total halo mass in a simulation with standard SF (KT07, left) and with H2-regulated SF 
(KMT09_ZFzlO, right) at z = 4. The dotted horizontal line indicates an equal to the cosmic baryon fraction fl\j/flM: i-e. a 100% gas to 
star conversion efficiency. 



modeling of stellar feedback processes in cosmological 
numerical simulations is beyond current computational 
capabilities, and it remains unclear how well the various 
sub-grid physics implementations discussed in the liter- 
ature actually capture the true nature of this feedback. 
Although we don't doubt that stellar feedback processes 
play an important role in regulating star formation in low 
mass halos, we now show that regulating SF by the H2 
abundance may play an equally important role in lower- 
ing stellar mass fractions in low mass (< 10^^ M©) halos, 
as previously suggested by G09, GKIO. 

The right panel of Figure [TT] shows that our H2- 
regulated SF prescription introduces a threshold below 
which the stellar mass content of halos is strongly sup- 
pressed. At halo masses greater than 10^^ M©, the star 
formation efficiency is only mildly reduced, dropping 
from (f^) = 0.062 to 0.035. Between 5 x 10^ and 
10^^ M0, only a small fraction of halos (about 10%) has 
been able to form stars and only with reduced SF effi- 
ciency, and at even lower halo masses star formation has 
been almost completely suppressed. We visually demon- 
strate this suppression in Figure [T2j where we show direct 
comparisons of the baryonic structure (S*, Sgas, ^Hs? 
and density- weighted metallicity) between the KT07 and 
KMT09_ZFzlO simulations for two representative halos; 
one high mass halo (Mh = 2.0 x 10^^ Mq) in which H2- 
regulation has not had a big effect, and one low mass halo 
(Mh = 5.8 X 10^ Mq), in which the stellar mass fraction 
has been suppressed by nearly one order of magnitude 
from f^ = 0.057 to 0.0069. Many more halos like these 
exist in our simulations. 

Figure [13] shows the baryonic content of halos in KT07 
and KMT09_ZFzlO as a function of their mass. We plot 
mean mass fractions in Mh-bins of width 0.25 dex of the 
total baryonic content, stars, neutral gas (HI + Hel + 
H2), ionized gas (HIT + Hell + Helll), H2, and metals. 
The total baryonic mass fraction remains equal to the 
cosmic mean down to ~ 2 x 10^ M©. The slight drop in 
the total baryonic and neutral gas fractions at lower halo 
masses can be attributed to the meta-galactic UV back- 



ground, which is able to ionize and heat most of the gas, 
preventing it from falling into the halo. The mass frac- 
tions of H2, metals, and stars roughly trace each other. In 
KT07, the stellar mass fraction decreases gradually from 
~ 0.4 of the cosmic baryon fraction at Mh = 10^^ Mq 
to 0.25 at 10^^ Mq, before the combined actions of UV 
background and (weak) supernova feedback further re- 
duce it to - 0.01 at 10^ Mq. In KMT09_ZFzlO, on the 
other hand, the stellar mass fraction is suppressed al- 
ready at much higher halo masses, dropping below 0.1 of 
the cosmic mean at 10^^ Mq and cutting off completely 
at 4 X 10^ Mq. Note that because we have included halos 
with M^ = in the calculation of the mean f^ it should 
not be viewed as a typical value for any individual halo, 
but as a population average. With H2-regulation the SF 
is quenched in low mass halos without heating and re- 
moving much of the gas, contrary to the effects of efficient 
supernova feedback. 

Why should a locally H2-regulated SF prescription be 
sensitive to the total halo mass? Recall that with the 
two-phase equilibrium assumption (KMT09 model), the 
local H2 abundance is completely determined by the HI 
column density and the metallicity of the gas. A reduced 
SF efficiency in low mass halos could thus be due to either 
lower column densities or to lower metallicities than in 
more massive halos, or to a combination of the two. 

Figure [14] reveals that lower metallicities are primarily 
responsible. The figure shows two-dimensional volumet- 
ric probability density functions (phase diagrams) of to- 
tal gas column density Ugas and metallicity Z, for / = 7 
grid cells in M< 10^° Mq (left panel) and M> 10^^ Mq 
halos (right panel) halos. The color scale indicates the 
volume fraction at a given (I]gas,Z), and the diagonal 
lines represent contours of constant /hs, according to 
the KMT09 prescription. The one-dimensional distri- 
butions at the top and sides of the figure show that 
while both low and high mass halos have comparable 
column density distributions, the metallicity distribution 
peaks at Z/Zq ^ 0.25 for the high mass halos, but at 
Z/Zq = 10-^ for the low mass halos. Recall that we 
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Figure 12. A comparison of the baryonic structure of two identical halos, one with high total mass (haloOOOl, left panel) and one with 
low total mass (haloOlSS, right panel), between the KT07 (left columns) and KMT09_ZFzlO (right columns) simulations a,t z = 4. From 
top to bottom, we show surface densities of stellar mass, total gas, and H2, and a density- weighted projection of metallicity, in a 5 x 5 kpc 
region centered on the galaxies. These are merely two representative halos, and many more like them exist in our simulations. The stellar 
content is greatly suppressed in the KMT09 low mass halo. 
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Figure 13. The baryonic content of halos versus their total mass, 
for the KT07 {top panel) and KMT09_ZFzlO (bottom panel) sim- 
ulations at 2; = 4. We show mean mass fractions normalized to 
the cosmic baryon fraction (Qi^/Qm) in M^-bins of width 0.25 dex 
for the total baryonic content (stars + all gas + metals; solid lines 
with circles), stars (dashed with circles), neutral gas (HI + Hel + 
H2; dot-dashed with upward triangles), ionized gas (HII + Hell + 
HeHI; dot-dashed with squares), H2 (dot-dashed with downward 
triangles), and metals (dot-dashed with diamonds). Note that al- 
though an H2 curve is plotted in the KT07 panel, H2 has no effect 
in that simulation. 



imposed a Z/Z© = 10~^ metallicity floor at z = 9, so 
this indicates that low mass halos typically remain un- 
enriched. This lack of metals in low mass halos is respon- 
sible for a reduced H2 abundance, and hence a suppressed 
stellar mass fraction. 

As we discuss in the remainder of this section, a num- 
ber of factors, both of physical nature and pertaining to 
the limitations of our simulations, affect the mass scale 
and the degree of the f^ suppression due to H2-regulated 
SF. A complete suppression of SF in halos with masses 
below 5 X 10^ M0 is too strong of an effect to reproduce 
the luminosity function of Local Group satellites, many 
of whi ch are consistent w ith having formed in < 10^ Mq 
halos (|Madau et al.|2008|). In this work, we demonstrate 
that H2-regulated SF can have important consequences 
for the star formation efficiency in low mass halos, but 
openly acknowledge that the true nature of this suppres- 
sion will depend on details of the star formation pro- 
cesses, stellar feedback, and metal enrichment that are 
not captured with sufficient fidelity in our simulations. 

4.1. Dependence on Two-Phase Equilibrium 



The assumption of two-phase equilibrium between 
a Cold Neutral Medium (CNM), hosting molecular 
clouds and SF, and the surrounding Warm Neutral 
M edium (WNM) is not uncontroversial. As discussed 



Krumholz et al. (2009), shocks from supernova explo- 



sions or supersonic ISM turbulence, for example, likely 
temporarily drive the surrounding gas out of pressure 
balance, in which case the typical gas density could dif- 
fer significantly from the value estimated from two-phase 
equilibrium. The time to re-establish equilibrium be- 
tween a CNM and WNM can be comparable (few Myr, 
Wolfire et al. (j2003|) to the shock recurrence time scale, 
and hence it is prudent to consider the effects of the 
KMT09 model without the two-phase equilibrium as- 
sumption. 

In this case the molecular hydrogen fraction is no 
longer independent of the Il2-dissociating LW back- 
ground, and absent a full radiation transfer treatment 
(G09) we need to externally specify its intensity. For 
computational ease we have considered here only the case 
of a constant, spatially uniform LW intensity. Of course 
it would be preferable to tie the LW intensity to the local 
SFR, averaged over ~ kpc scales and tens of Myr, but we 
defer the investigation of this more realistic treatment to 
future work. 

In Figure 15] we show (f^) vs. Mh at z = 5 for simu- 
lations KMTU9_FLW{1, 10,100,1000}, with values of the 
LW intensity ranging from 1 to 1000 times the z = 
Milky Way value of 7.5 x 10~^ LW photons cm~^ (Draine 
^1978} . For clarity we have plotted the mean values of 
in bins of Mh of width 0.25 dex. The means include 
halos with = 0, which implies that these values of 
(f^) should thus be viewed as population averages, rather 
than representative values for any individual halo. The 
results from the standard SF (KT07) and the two-phase 
equilibrium (KMT09) simulations are included for com- 
parison. The lowest LW intensity case (KMT09_FLW1) 
is almost indistinguishable from the standard SF model. 
In this case the dissociating flux is not strong enough 
to affect the II2 abundance even in the lowest mass ha- 
los. As the intensity of the LW background is increased, 
the halo mass at which H2 formation begins to be sup- 
pressed increases in proportion. At 1000 times the Milky 
Way's LW intensity (lowest dotted line), the suppression 
mass almost reaches the value in the two-phase equilib- 
rium model. Recall that the KMT09_FLW simulations 
were run with an effective sub-grid clumping factor of 30 
(((Jc^^_2i/7^-i6.5) = 1/30), and the resulting enhanced 
H2 recombination rate is partially responsible for the re- 
duced SF suppression compared to the two phase equilib- 
rium simulation. However, while the absolute strength of 
the SF suppression in the simulations without two-phase 
equilibrium is sensitive to the value of the clumping fac- 
tor, the trend with the intensity of the LW background 
is not. 

4.2. Resolution Dependence 

In Figure 16 we compare f^ vs. halo mass for the 
KMT09_L8 simulation to the results from KT07 and 
KMT09. While there still is a suppression of star for- 
mation in low mass halos compared to KT07, the scale 
at which this suppression occurs has shifted to lower 
masse s, ar ound 3 x 10^ Mq. Considering our earlier find- 
ing (§|3.4[) that the additional resolution increases the SF 
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Figure 14. 2D phase diagrams of Sgas vs. Z determined at /max = 7 in the KMT09 simulation at z = 5, for low mass (M< lO-*^^ M©, left 
panel) and high mass (M> 10^^ Mq, right panel) halos. Contours of constant (at 0.1, 0.5, and 0.9) are plotted with solid lines. For 
reference we show with dotted lines the column densities corresponding to the 3D density thresholds employed in the KT07 (50 cm~^), 
KT07_low (5 cm~^), and KT07_high (500 cm"*^) simulations. The majority of gas in low mass halos has very little or no H2, and this 
appears to be primarily due to lower metallicities, not lower column densities. 
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Figure 15. {f^> vs. Mh at z = 5 for the four KMT09 simulations 
without the two phase equilibrium assumption and a spatially uni- 
form Lyman- Werner radiation field with intensity equal to 1, 10, 
100, 1000 X the mean Milky Way value (dotted lines with down- 
ward triangles; increasing LW intensity from top to bottom). We 
plot the mean of f^ in M/^-bins of width 0.25 dex. The means 
include halos with f^ = 0. For comparison we also show the orig- 
inal KT07 (solid line with red circles) and two-phase equilibrium 
KMT09 results (dashed line with upward triangles). Note that 
for the KMT09_FLW simulations (but not for KMT09) we apply 
a "clumping factor" of 30 to the H2 formation rate (see text for 
details). 



rates in KMT09_L8, it is perhaps not surprising that the 
f^-suppression mass scale also exhibits some resolution 
dependence. In this case the sensitivity to the maximum 
resolution arises because the increase in density afforded 
by the additional refinement level exceeds the factor of 
two reduction in grid cell width, thus leading to a net 
increase in column densities on the finest grid cells. The 
enhanced LW shielding eases the transition to the molec- 
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Figure 16. f^ vs. halo mass (at z = 6) for the KMT09_L8 sim- 
ulation, which has one additional refinement level compared to 
KMT09 and KT07. The halo mass scale at which star formation is 
suppressed due to the inability of gas to become molecular is lower 
than in the KMT09 simulation. The true value of this suppression 
mass scale in nature will depend on the strength of stellar feedback, 
metal mixing, and other processes not resolved in our simulations. 



ular phase, and allows star formation to occur in lower 
mass halos. 

The exact value of this transition will depend on the 
details of the star formation and stellar feedback pro- 
cesses on scale s be low our current resolution limit, as we 
discussed in § |3.4| However, our fiducial /max = 7 grid 
resolution of 76.3 x 5/(1 + proper parsec is compara- 
ble to the observed size of giant atomic-molecular cloud 
complexes, and our / = 7 grid cell averaged densities 
left panel) are in good agreement withobser- 



(Fig. 10 



vat ionaf estimates of their average densities (Blit z|1993 ). 
We thus believe that the resolution of our fiducial simu- 
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Figure 17. The dependence of (f^(Mh)} at z = 6 on the metal- 
Ucity floor. We plot the mean of in M^-bins of width 0.25 
dex at z = 6 for KT07 (sohd hnes with circles) and three ver- 
sions of KMT09 with increasing values of the z = 9 metallicity 
floor, [Zfloor] = —3.0 (dashed with upward triangles), —2.5 (dot- 
ted with downward triangles), and —2.0 (dot-dashed with squares), 
and one case with [Zfloor] = —3.0 applied at z = 10 (dashed with 
diamonds). The means include halos with = 0. Note that a 
simulation with log2o(Zfloor/Z0) = —4.0 at 2; = 9 produced no 
stars at all. 



lations 
studying. 



7) is well matched to the problem we are 



4.3. Metallicity Floor Dependence 

In Figure [Tt] we show how the suppression depends 
on the amplitude and time of the metallicity floor that we 
impose early in the simulations to mimic the enrichment 
from Pop. Ill supernovae. Higher values of [Zfloor] and 
an earlier time allow SF to occur in lower mass halos. 
Note that the KMT09_ZF4.0 simulation ([Zfloor] = -4.0) 
produced no stars at all. Imposing the metallicity floor at 
z = 10 instead of z = 9 has the particularly interesting 
result of allowing a small amount of star formation in 
very low mass halos (Mh < 5 x 10^ Mq), which may help 
to alleviate the tension with the faint-end of the Local 
Group satellite galaxy luminosity function. 

For a given Zfloor, a halo's ability to form stars is de- 
termined by the highest column densities that its gas can 
condense to prior to the enrichment from the first star 
particles in our simulation. Gas in more massive halos 
is able to reach higher column densities, owing to their 
earlier collapse times and deeper potentials. Below some 
halo mass scale, gas will not be able to reach sufficiently 
high column densities to turn molecular and allow star 
formation to occur. The necessary column densities de- 
pend on metallicity: at [Z] = log2o(Z/Z0) = —3.0, gas 
must reach a column density of 6.7 x 10^ (1.4 x 10^, 
8.0 X 10"^) M0 pc-2 in order to be 10 (50, 90) per cent 
molecular. At [Z] = —2.5 the required column is reduced 
to 2.3 X 10^ (4.9 X 10^, 3.0 x 10"^) M© pc-^ and it is only 
8.0 X 10^ (1.7 X 10^ 1.0 X 10"^) M0 pc-2 at [Z] = -2.0. 
This naturally explains the observed dependence of the 
f^ suppression scale on the amplitude of the metallicity 
floor. It implies that the true f^ suppression realized in 
nature will depend on the details of the enrichment his- 
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Figure 18. The luminosity function in the KT07 (open circles) 
and KMT09_ZFzlO (solid triangles) simulations compared to obser- 
vational results. We calculate UV luminosities from the simulated 
galaxies' SFR according to Luv = 8.0 x 10^'^ (SFR/M© yr-^) erg 



Hz- 



(Mt 



51.63 — 2.5 log ]^Q(Luv/6rg s ^ Hz ^)), c orre- 



sponding to a Salpeter IMF from 6.1 to 125 Mq (same as Bou wens| 
let al.'2011'). Error bars are statistical only. The solid lines are the 
uncorrected luminosity functions reported by Bouwens et al. ( 2 007| 
2011 ), and the dashed lines the same relations corrected for dust ex- 
tmction by fl.5 5, 0.625, 0.375, 0, 0) m agnitudes at z = (4, 5, 6, 7, 8) 
(cf. Table 8 in [Bouwens et al.|201lt . 



tory of a given halo, which should lead to a broadening 
of the Mh dependence of the f^ suppression. 

5. COMPARISONS WITH HIGH REDSHIFT OBSERVATIONS 
5.1. Evolution of the Luminosity Function 



In Figure 18 we present a comparison of the luminos- 
ity functions(LF) from the KT07 and KMT09_ZFzlO 
simul ations and the recent hi^h re dshift determinations 



from Bouwens et al. (2007 2011) based on deep HST 



ACS and WhV'i observalions oTS {{z) = 3.8), V (5.0), i 
(5. 9), z (6.8), and F - band (8.0) dropout galaxies. Follow- 
ing |Bo^ens~er^D (|20^ hereafter BIO), we determine 
UV luminosities trom the simulated SFR according to 
Luv = 8.0 X 10^^ (SFR / M© yr-^) erg s"^ Hz-\ cor- 
responding to a Salpeter IMF fr om 0.1 to 125 M^;^ and a 
constant SFR over > 100 Myr (|Madau et aT] |1998D . We 
calculate SFR from our simulated galaxies by summing 
the mass of all young star particles with ages less than 
= 30 Myr and dividing by the SF time scale. 



SFR 



E 



(11) 
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The assumption of a constant SFR over the past 100 Myr 
is not hkely to hold for our simulated galaxies, and we 
therefore may be overestimating their UV luminosities by 
a factor ^ 2 or so. A Kroupa IMF, on the other hand, 
would result in ~ 1.7 higher UV luminosity for a given 
SFR. We compare the simulated LFs to the Schechter 
function fits reported by BIO (shown in solid lines), as 
well as to these same relations corrected for dust extinc- 
tion (dashed lines) by (1.55, 0.625, 0.375, 0, 0) magni- 
tudes at z = (4, 5, 6, 7, 8) (cf. Table 8 in BIO). 

The KT07 LF exceeds the uncorrected BIO LF by - 1 
dex at all redshifts less than 8, and over the entire range 
of luminosities probed by BIO. Applying the BIO dust 
correction brings the z = 4 LF into agreement with the 
KT07 simulation. This agreement is remarkable, since 
it results without any tuning of our models. It is puz- 
zling, however, since it seems to imply that there is no 
dwarf galaxy problem ai z = 4. In general, the steep 
faint-end LF slopes (—1.7 to —2.0) at z > 4 reported by 
BIO appear to be in tension with the need to strongly 
suppressed star formation efficiency in low mass halos 
in order to match the lower redshift stellar mass func- 
tions and dwarf galaxy luminosity functions. Regardless 
of this, at z > 4 the disagreement between our KT07 
LFs and the dust-corrected observational ones remains 
substantial. 

The KMT09_ZFzlO LF matches the KT07 one at high 
luminosities, reflecting the fact that H2-regulation is not 
effective in high mass (high Luv) halos. Close to the sen- 
sitivity limit of the BIO observations, the KMT09_ZFzlO 
LF begins to roll over, improving the agreement between 
simulations and observations. This rollover, however, 
continues to lower luminosities, and it appears that in 
the KMT09_ZFzlO simulation the H2-suppression may 
be too efficient at low UV luminosities, at least if the 
faint end of the observed UV LF continues to rise to 
lower luminosities. As we discussed in § 4.2, the exact 
mass scale of the SF suppression (and hence the down- 
turn of the LF) is sensitive to the nature of the metal 
enrichment process, and also depends on the resolution 
of our simulations. 

5.2. Evolution of the Stellar Mass Density 



In Figure 19 we show the redshift evolution of the co- 
moving stellar mass density (SMD) contributed by galax- 
ies with > lO^ '^^ Mo in the KT07 and KMT09_ZFzlO 
simulations, compare d to t he ob servational determina- 
tion from Gonzalez et al. (2011). These authors used 
restframe optical photometry from Spitzer/IRAC to infer 
stellar masses for a sample of ~ 300 z = 4 galaxies, and 
combined these with HST ACS and WFC3/IR restframe 
UV fluxes to establish an empirical stellar mass to UV lu- 
minosity relation. Assuming this same M^-Luv extends 
to higher redshifts, they integrated the |Bouwens et al.| 
( [2011^ z = 4 — 7 restframe UV luminosity function down 
to a limiting magnitude of Muv = — 18 (corresponding 
to M^ = 10^"^^ M0), and obtained an empirical estimate 
of the SMD evolution: SMD ex (1 + ^)-3-4±o.8^ shown as 
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To demonstrate the systematic uncertainty in this kind 
of empirical determination of the stellar d ensity evolu- 
tion, we also show results from [Labbe et al. (|2010a|) (dot- 
ted line, SMD oc (1 + z)-^) and ILabbe et all (1201 Ob) 
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Figure 19. The stellar mass density for galaxies with > 
10'^ '^^ M0 vs. z in the KT07 and KMT09_ZFzlO simulations 
compared to observational determinations from [Go nzalez e t al.| 
(2011) (SMD oc (l + z)-3-4±0-8, solid line), [Labblet al. {20fn^ 
(SMD oc (1 + z)-^, dotted line), and Labbe et al. (2010b) (SMD 
oc (1 + 2;)~^-^=*='-'-^, dashed line). The uncertaintv i n the s lope re- 
ported by Gonzalez et al. (2011) and Labbe et ah] ( [2010b| is indi- 
cated by the gray regions, arbitrarily normalized~at z = 4. 

(dashed line, SMD cx (1 + 2;)~^"^^^"^), who performed 
similar work including z = S Y-band dropouts, but used 
M^-Lxj v relations derived from different low redshift data 
dStark et al. 2009). 

Likewise our simulated SMD's are not free from sys- 
tematic uncertainties: on the one hand they are likely 
somewhat underestimated, especially at high redshift, 
since the simulations' limited box size (12.5 Mpc) re- 
sults in a delayed formation of massive halos, owing to 
the absence of density perturbations in the initi al condi- 
tions wit h wavelengths exceeding the box size (Tormen 
fc Bertschinger |l996 ). On the other hand the lack of 
ettective stellar feedback in our simulations probably ar- 
tificially enhances the stellar mass density at all redshifts. 

Given the substantial uncertainties in both observa- 
tions and simulations, we consider the agreement be- 
tween the two to be satisfactory. We do note, however, 
that the simulati ons prefer a somewhat steeper evolution 
of the SMD than |Gonzalez et al.| ( [2QTT| , and that the H2- 
regulated SF reduces the SMD by a factor ~ 3— 5 over the 
covered redshift range compared to the standard SF case. 
The difference between KT07 and KMT09_ZFzlO de- 
creases towards lower redshifts, which can be attributed 
to the increase of the typical halo mass (the knee in the 
Press- Schechter massfunction) with decreasing redshift, 
resulting in an increasing fraction of resolved halos with 
mass above the f^-suppression mass scale. 

5.3. Evolution of the Star Formation Rate Density 

In Figure |20] we extend the comparison between our 
simulations and high-z observational data to the evolu- 
tion of the SFR density, obtained by summing the SFR 
(Eq. [n]) of ah galaxies with M^ > lO^ '^^ M© (Mw < 
— 18) and dividing by the box volume. The observational 
data again comes from BIO, who derived SFR densities 
by converti ng their restframe U V luminosities to SFR ac- 
cording to Madau et al. (1998), and integrated down to a 
limiting magnitude of Muv = —18, with (orange band) 
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Figure 20. The SFR density for galaxies with > W-^^ M© 
vs. z in the KT07 and KMT09_ZFzlO simulations, compared to 
observational determinations from Bouwens et al. (2011), which 
were obta ined by converting restframe U V luminosities to SFR 
(|Madau et al. 1998), without (blue band) and with dust corrections 
(oran ge band), and by conve rtin g stellar mass densities (Stark et al. 
|2009 Labbe et al.|2dl0a|b[ [Gonzalez et al.|201l| > to Sl^'H densities' 
(light gray bancT^ 

and without (b l ue ban d) a dust correction at z < 6 from 
Bouwens et al. (|2QQ9 ). Note that the ver y blue UV con- 
tinuum sl opes observe d in z > 7 galaxies (Bouwens et al. 



25TU1 [Bmiker et al. 2QTq1 [Fmkelstein et al. 2U10 ; OescE 
"^^|2Qir' ' 



et al. imply low dust abundances and correspond- 

ingly weak extinction corrections at these high redshifts. 
For comparison we also show (light gray band) the SFR 
density implied by the SMD evolution discussed above, 
assuming a fixed SMD cx (1 + z)~^-^ extrapolation to 
z>S (Labbe et al. 2010b). 

Our simulated S FR densities are s omew hat higher than 
those reported by Bouwens et al. (2011). H2-regulated 
SF, however, reduces the Sl^'R densities by about a factor 
of two, and brings them close to agreement with those 
determined from the SMD evolution, and only a factor 
~ 3 above the luminosity density derived values at 2; > 6. 
At lower redshifts, the rise due to dust corrections in the 
observational SFR densities reduces the difference, and 
an extrapolation of the KMT09_ZFzlO SFR densities to 
z = 4 is in good agreement with the dust corrected obser- 
vational determination. If the observed steep UV contin- 
uum slopes are a selection effect and not representative of 
typical z >7 galaxies, then a more heavily dust obscured 
population of galaxies that has been missed by current 
surveys could close the high redshift gap between our 
simulations and current data. Additionally it is possible 
that a more effective stellar feedback implementation in 
our simulations could further reduce the simulated SFR 
densities. 

6. SUMMARY AND DISCUSSION 

Motivated by the observational fact that star forma- 
tion, both in the local universe (e. ^. Bi^el et al . 2008 j) 
and at intermediate redshifts (e.g. [Genzel et al. |5oiO ), 
correlates more tightly with the density ot gas in the 
molecular phase than the total gas density, we have im- 
plemented an H2-regulated SF prescription in cosmolog- 
ical galaxy formation simulations, with the goal of inves- 



tigating whether this new ingredient can help to alleviate 
outstanding problems in our theoretical understanding of 
dwarf galaxy formation. 

Rather than following the non-equilibrium molecular 
hydrogen chemistry including time-dependent and spa- 
tially inhomogeneous radiation transfer of the ionizing 
and H2-dissociating stellar radiation field (G09), we uti- 
lize in our simulations the results of one-dimensional 
radiative transfer calculations of the H2 formation- 
dissociation balance in an idealized spherical giant 
atomic-molecular complex subject to a uniform and 
isotropic Lyman- Werner (LW) radiation field ( |Krumholz| 

] [2QTq1 ). I'hese cal- 



et al. 



2008 



2009 McKee & Krumholz 



culations showed that under the assumption of two-phase 
equilibrium in the ISM, the H2 abundance is determined 
entirely by the column density and metallicity of gas on 
~ 100 pc scales - comparable to the cell size of our high- 
est resolution grid cells. 

We have deliberately chosen a weak implementation 
of supernova feedback in order to isolate the effects of 
the molecular hydrogen chemistry on star formation in a 
cosmological setting. Remarkably, we find that account- 
ing for the transition from HI to H2 alone is able to fulfill 
one of the functions of the neglected supernova feedback, 
namely the suppression of star formation in low mass ha- 
los. The main results of our study are summarized as 
follows. 

1. Both our conventional and H2-regulated SF pre- 
scriptions are able to reproduce the observational 
scaling relation between the SFR surface density 
and the total gas density, the Kennicutt- Schmidt 
relation (Fig. [3|. 

2. With the conventional SF prescription the observed 
turnover of Ssfr at low Sgas is recovered by man- 
ually tuning a SF threshold density. The H2- 
regulated prescription, however, automatically re- 
produces this cutoff (Fig. |5| and thus reduces the 
number of free parameters in the SF prescription 
by one. 

3. We are able to reproduce many of the observa- 
tional results pertaining to the KS relation, for 
example the slope of the molecular gas KS rela- 
tion (SsFR — ^Hs) and the much weaker correla- 
tion of SsFR with the atomic gas surface density 
(Fig. [6|. We also recover the observed metallic- 
ity dependence of the low Sgas turnover in the KS 
relation (Fig. which occurs at higher Sgas in 
lower metalliciiy systems, refiecting the atomic-to- 
molecular gas transition. Lastly, in agreement with 
recent observations, we see an increased scatter in 
the H2 KS relation for smaller spatial smoothing 
scales. 

4. We expect an increased scatter in the KS relation 
at higher redshifts, since the typically larger stellar 
velocity dispersions in high redshift systems will 
allow young stars to wander outside of high column 
density regions. 

5. H2-regulation suppresses SF in low mass halos, re- 
ducing the need for stellar feedback (Fig. 11), as 



previously reported by G09 and GKIO for a set of 
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cosmological zoom-in galaxy formation simulations 
including non-equilibrium H2 formation and radia- 
tive transfer. We confirm these prior results with 
better statistics (hundreds of halos) and without 
the need for a complicated and expensive radiative 
transfer treatment. 

6. The halo mass dependence appears to be tied pri- 
marily to a difference in metal enrichment, rather 
than gas column density (Fig, [m] ). On the other 
hand, for a given metallicity floor (set for exam- 
ple by the first generation of Pop. Ill supernovae), 
the star forming ability of halos is determined by 
the highest column densities their gas can condense 
to. Low mass halos don't have sufficiently high col- 
umn density gas to allow the transition to molecu- 
lar phase and hence star formation. 

7. Suppressing the star formation efficiency in low 
mass halos lowers the cosmic stellar mass density 
and star formation rate density. We find reason- 
able agreement between our H2 -regulated simula- 
tion and observational determinations of the evo- 
lution of the SMD and SFR-density at z > 4. Both 
simulations and observations are subject to large 
systematic uncertainties. 

As we discussed throughout the text, a number of 
caveats apply to our findings. First, the results we have 
presented here are not completely independent of the 
implementation details of our simulations (e.g. 2-phase 
equilibrium assumption, metallicity floor, numerical res- 
olution). We therefore cannot claim with any confidence 
to have pinpointed the halo mass below which star for- 
mation is suppressed due to the inability of gas to become 
molecular. 

Second, given our computational resources, we were 
unable to continue our simulations beyond z = 4. It is 
possible that a large build-up of atomic gas will lead to 
a burst of star formation at lower redshift {z o:^ 1 — 3), 
when metallicities have increased sufficiently to allow the 
transition to H2 to occur, and in fact this may help to 
explain how the specific SFR of 10^° - 10^^ M© halos 
can excee d the instantaneous gas a ccretion rate at these 
redshifts (iKrumholz fc Dekel|[2QTT| . 

Lastly, like most cosmological galaxy formation simu- 
lations to date, our simulated galaxies suffer from the so- 
called "baryonic overcooling problem", resulting in un- 
realistically high central densities (and hence strongly 
peaked circular velocity curves) and stellar mass frac- 
tions in our high mass halos that are too large compared 
to observations. Although we cannot exclude the possi- 
bility that these artificially high central densities affect 
the results of our study, we believe that this makes our 
results conservative, in the sense that a more realistic 
simulation with lower central gas densities would likely 
experience even more Il2-regulated SF suppression. Re- 
duced central densities would demand higher metallicity 
for the transition to II2 to occur, and at the same time 
they would allow the metals ejected by stars to become 
diluted more easily, thereby reducing the metallicity of 
the high density star- forming gas. 

Even allowing for these caveats, we believe that the 
atomic-to-molecular gas transition may play an impor- 
tant role in regulating star formation in low mass halos. 



It may also help to expla i n the result recently reported 
by Boylan-Kolchin et al. (2011) that collisionless simu- 
lations of galactic dark niatter substructure (Via Lactea 
II and Aquarius) predict a number of subhalos too cen- 
trally concentrated to host any of the kn own Milk y Way 
dwarf satellite galaxies. As Boylan-Kolchin et al. (2011) 
suggest, one possible explanation of this puzzle is that 
star formation in dark matter halos becomes stochastic 
below some mass. The sensitivity of an H2-regulated SF 
prescription to early metal enrichment may provide the 
necessary stochasticity. 

In our simulations, H2-regulation quenches SF in dwarf 
galaxies at the outset, without the need to artificially en- 
hance stellar feedback by turning off gas cooling or hy- 
drodynamically decoupling momentum-driven winds, as 
is commonly done in the literature. These latter imple- 
mentations of enhanced stellar feedback quench SF in 
dwarf galaxies due to the inability of their shallow po- 
tential wells to retain the resulting gas outflows. In H2- 
regulated SF, by contrast, it is the inability of low mass 
halos to accumulate a sufficiently large column density 
of metal-enriched material that results in suppressed SF. 
Efficient supernova feedback implies the presence of a 
large pool of hot gas as well as galactic winds, neither of 
which are necessarily present in an H2-regulated scenario. 
This difference may be be a good way to observationally 
distinguish the two mechanisms. 

Taken at face value, our results imply that many low 
mass dark matter halos at high redshift should be filled 
with relatively cold, yet atomic gas, which is prevented 
from becoming molecular by its low metallicity. While 
neutral gas surveys with radio telescopes are hard pressed 
to reach beyond the local universe, constraints on the HI 
mass funct ion do exist at low redshift. The ALFALFA 
survey (Martin et al.ll2010l), for example, has reported a 
z < 0.06 HI mass function down to log(MHi/M0) = 6.2, 
that is well fit by a Schechter function with a low mass 
slope of a = —1.33. This is considerably shallower 
than the mass func tion of dark matter halos (a ~ —1.8, 
|Boylan-Kolchin et al. 2009), indicating that not all low 
mass dark matter halos can be allowed to retain their 
full baryonic content in the form of cold, atomic gas 
down to the present epoch. Another constraint may come 
from the s tatistics of damped Lyman-a (DLA) systems 



at z < 5 (Prochaska fc Wolfe 
20091). Whether the observed 



2009 Noterdaeme et al. 



requency distribution of 
DLA column densities and their covering fraction is con- 
sistent with large amounts of atomic gas in high redshift 
dwarf galaxies should be investigated in detail. 

Strongly suppressing star formation in low mass halos 
at high redshift could make it difficult for faint galaxies to 
reionize the univer se by z = 6, as is commonly advocated 
(e.g. Madau et al. 1999). Indeed, in our simulations the 
total number of hydrogen ionizing photons per hydrogen 
aton|^ produced by z = 6 is reduced from ^ 18 in KT07 
to ^ 6 in KMT09. As shown by rtKuhlen fc Faucher 



[Ciguere |2012|), a suppression scale of Mh = 10^^ M© 
can simultaneously satisfy reionization and lower red- 
shift Lyman-a forest constraints only if the escape frac- 

^ We assume 4000 hydrogen ionizing photons per stehar baryon 
for every solar mass of stars formed, as appropriate for a Salpeter 
IMF from 0.1 to 125 and a mean ionizing photon energy of 20 
eV ( [Madau et al.|1999[ [Leitherer et al.|1999| >. 



21 



tion of ionizing radiation evolves strongly from z = 4 to 
higher redshifts. A slightly smaller supp ression mass of 

10^ M0 (e.g. as in'Kru mholz fc Dekel|[2Qll| ), however, 
is in good agreement with all constraints. 

While we have emphasized in this work that H2- 
regulated SF can perform some of the functions typi- 
cally assigned to supernova feedback, we of course ac- 
knowledge that supernovae do in fact occur in nature 
and that their associated injection of energy and momen- 
tum into the surrounding interstellar medium is likely to 
significantly impact subsequent star formation, and may 
help to solve the problem of forming galaxies that are 
too concentrated in the simulations. Future research will 
be necessary to elucidate how the interplay of molecular 
chemistry and supernova feedback shapes star formation 
in dwarf galaxies and beyond. 
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APPENDIX 
SOBOLEV-LIKE APPROXIMATION OF E 
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Figure 21. Left: Comparison of a Sobolev-like estimate of the column density (X^sob = p x (p/Vp), see text for details) with the simple 
grid-cell column {^ce\\ = pAx), obtained by post-processing the z = 4 KMT09 output. Only columns greater than 50 Mq pc~^ are shown, 
since the H2 content is negligible at lower columns. The solid and dashed black lines show the median and 16*^ - 84*^ percentiles of Sgob 
for a given Sceii- The Sobolev-like estimate results in slightly lower columns, especially at high columns. Right: KS relation for simulation 
KMT09_Sob, in which the Sobolev-like approximation was used at simulation run-time to calculate the column density entering the KMT09 
expression for fn^. Note that we did not use the Sobolev-like approximation to determine the x-axis quantity Sgas, the column smoothed 
at ^ = 3 (~ 1 kpc). The slightly lower columns with the Sobolev-like approximation result in lower fn^ and hence a slightly reduced SPR 
at a given Egas- 



One way to alleviate the explicit resolution dependence in our implementation of the KMT09 algorithm (Eq. [4j[9|) is 
to use a Sobolev-like approximation of the column density (see Gnedin et al. 2009) , instead of simply multiplying the 
grid cell density by its width. We have implemented a Sobolev-like approximation as follows: 
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and similarly for the y and z directions. This definition has the virtue that it is resolution-independent, at least to the 
extent that the density field itself is resolution-independent. To obtain a total column density, we take the harmonic 
mean over the six cardinal directions. 



1/E+ + + 1/E+ + 1/E^ + 1/E+ + 1/E^ 



(A2) 



We have applied this Sobolev-like approximation for Sgas in a post-processing analysis of the z = 4 output of the 
KMT09 simulation, and compared the resulting column density with the simple grid cell based estimate used in the 
simulation, see the left panel of Fig. 
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We're only showing points with S > 50 Mq pc~ , since at lower columns the 
H2 abundance, and hence SFR, is negligible. The median of Ssob lies close to the 1-1 line, but drops below it around 
Sceii = 100 M0 pc~^, implying that the Sobolev-like approximation systematically yields slightly lower values at high 
Sceii. At Sceii = 100, 1000, and 5000 M© pc-^ the median of Ssob is 101, 880, and 3930 M© pc-^ respectively. The 
1-(T scatter of Ssob around the median is 0.26 dex (i.e. 68% of all cells lie within a factor of ~ 1.8 of the median). 
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Figure 22. Comparison of the KS relation for the KMT09_Sob and KMT09_SobL8 simulations at z = 6. The use of Esob does not fully 
remove the resolution dependence, but reduces it mildly compared to simulations with ^cell (cf- right panel of Fig. |10| ). 

This post-processing analysis indicates that a simulation actually using the Sobolev-like approximation at run-time 
might have somewhat reduced SF, owing to the systematically slightly lower columns, which provide less shielding 
and hence lower H2 abundances. This is in fact borne out in practice, as shown in the right panel of Fig. [2l] We 
ran two additional simulations, KMT09_Sob and KMT09_SobL8, in which the column densities entering the KMT09 
/hs prescription were determined at run-time using the Sobolev-like estimator. Indeed the resulting KS relation has 
a slightly lower normalization, indicating a somewhat reduced SF efficiency. 

As discussed in the main text, to some degree there will always be a resolution dependence in the determination of 
the column density in our simulations, simply because we are not resolving the true Jeans length and are applying a 
minimum pressure support in order to prevent artificial fragmentation. This is demonstrated in Fig. [22j which shows 
the KS relations for the KMT09_Sob and KMT09_SobL8 simulations. The corresponding plot for simulations with 
Hceii is the right panel of Fig. [Toj and a comparison between these two figures shows that the use of Ssob only mildly 
reduces the resolution dependence. 

The two different methods of estimating column densities (Sceib ^Sob) are resolution-dependent in different ways. 
The difference between the two in S and in the KS relation is fairly small (less than a factor of 2), and does not impact 
the overall conclusion of our study that the metallicity-dependent nature of the atomic-to-molecular transition can 
play a major role in explaining the low star formation efficiency in dwarf galaxies. 



